Nonthermal antiferromagnetic order and nonequilibrium criticality in the Hubbard model 



O 
o 

O: 



i: 



I 

O 

o 



> 

cn 

o 
d 



X: 



Naoto Tsuji,' Martin Eckstein,^ and Philipp Werner' 

'Department of Physics, University of Fribourg, 1700 Fribourg, Switzerland 
^Max Planck Research Department for Structural Dynamics, 
University of Hamburg-CFEL, 22607 Hamburg, Germany 
(Dated: October 5, 2012) 

We study dynamical phase transitions from antiferromagnetic to paramagnetic states driven by an interaction 
quench in the fermionic Hubbard model using the nonequilibrium dynamical mean-field theory. We identify two 
dynamical transition points where the relaxation behavior qualitatively changes: one corresponds to the thermal 
phase transition at which the order parameter decays critically slowly in a power law oc f^''", and the other is 
connected to the existence of nonthermal antiferromagnetic order in systems with effective temperature above 
the thermal critical temperature. The frequency of the amplitude mode extrapolates to zero as one approaches the 
nonthermal (quasi-)critical point, and thermalization is significantly delayed by the trapping in the nonthermal 
state. A slow relaxation of the nonthermal order is followed by a faster thermalization process. 
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In many physical systems out of equilibrium, phase transi- 
tions occur as a real-time process of symmetry breaking and 
symmetry recovery. Examples for such "dynamical phase 
transitions" include the evolution of the universe [1], liquid 
helium and photoinduced phase transition in solids [3- 
[st]. The macroscopic aspects are often described by the time- 
dependent Ginzburg-Landau theory, where the order parame- 
ter is supposed to vary sufficiently slowly in time and space, 
so that the system can be considered to be locally close to 
thermal equilibrium. On the other hand, recent experimen- 
tal developments of time-resolved measurement techniques in 
solids [6] and cold atoms ^ allow one to study dynamical 
phase transitions very far from equilibrium on the microscopic 
time scale of correlated quantum systems. In these cases, a 
'near-equilibrium' description might not be applicable. For 
instance, it has been recently suggested that superconductiv- 
ity can be induced above the equilibrium critical temperature 
(Tc) by coherently exciting certain lattice vibrations, and that 
it lasts for a relatively long time (a few tens of ps) before 
thermalization occurs |5]. This observation is reminiscent of 
the prethermalization phenomenon llil- iTlll . or the dynamics in 
the presence of a nonthermal fixed point in relativistic quan- 
tum field theories lll2ll . A fundamental question that we pose 
here is if the existence of such a nonthermal fixed point in 
correlated condensed matter systems allows symmetry broken 
states to survive above Tc, and how it affects the dynamics. 

An important and still unresolved issue is how to charac- 
terize a nonequilibrium phase transition and its critical be- 
havior for quantum systems |13, 14]. Previous studies have 
in particular focused on the dynamics near quantum pha se 
transitions in low dimensional systems (e. g., Refs. fllsMlSh . 
Higher dimensional systems are usually expected to show a 
thermal criticality out of equilibrium since quantum fluctua- 
tions are well suppressed. In this paper, we study a dynam- 
ical phase transition for a simple microscopic model of cor- 
related fermions, namely the Hubbard model. In equilibrium, 
the model exhibits a phase transition from paramagnetic (PM) 
to antiferromagnetic (AFM) order [see the phase diagram in 



Fig.lTJa)]. By changing the interaction in time, we cross the 
phase boundary dynamically. In particular, we explore the 
weak-coupling regime of the Hubbard model (fo r the strong- 
coupling side, see our complementary work [20]). Contrary 
to the naive expectation, we find that the nonequilibrium re- 
laxation behavior can be very dififerent from the thermal one 
even in the large-dimensional limit. A new phenomenon that 
we demonstrate here is that in addition to the thermal criti- 
cal point there exists one more quasi-critical point (or sharp 
crossover) at which some time (energy) scale almost diverges 
(vanishes). Between these points, the system is trapped in a 
nonthermal "ordered" state [Fig. [Ub)], where the order pa- 
rameter stays nonzero even though the effective temperature 
(which will be defined below) is above Tc- 
The model Hamiltonian is given by 
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where is the band dispersion, (c^o-) is a creation (an- 
nihilation) operator of fermions with spin <j, U is the (time- 
dependent) interaction strength, and n,o- = cj^c,o-. For con- 
venience, we take a semicircular density of states, D(e) = 
v/4 - (e/f*)2/(27rf*), and use t* (f*"') as the unit of energy 
(time). We only consider the half-filling case. The initial 
state is in thermal equilibrium with temperature T, which is 
chosen such that the initial value of the staggered magnetiza- 
tion m = ( |n-t- — n^l ) is 0.4. The interaction is changed as 
U(t) - Ui + (Uf - Ui)t/tg (0 < f < t^) with quench time 
t^j = 8 fixed. The interaction quench can be implemented in 
cold atom systems with the use of the Feshbach resonance, or 
by modifying the depth of the lattice potential, and has also 
been proposed to be possible in solids driven by strong elec- 
tric fields 1 21, 22]. 

The time evolution of the Hubbard model with AFM or- 
der is studied with the nonequilibrium dynamical mean-field 
theory (DMFT) ['23, 24]. It becomes exact in the large dimen- 
sional limit [25], where the self-energy becomes local in space 
but keeps dynamical coiTelations. When one allows for AFM 
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FIG. 2: (Color online). Time evolution of in for quenches (a) 
Ui = 2 ^ Uf = 1.0, 1.1,..., 1.9 (from bottom to top), and (b) 
Ui = 2.5 — > Uf = 1.5, 1.6, ... , 2.4. The arrows indicate the corre- 
sponding thermal values mth reached in the long-time limit. 



FIG. 1 : (Color online), (a) Equilibrium phase diagram of the Hub- 
bard model in the weak-coupling regime at half-filling, calculated by 
DMFT with several different impurity solvers. QMC data are taken 
from Ref. 19. Effective temperatures for quenches from a fixed initial 
state (Ui = 2, black dot) to various final states (open dots) are shown. 
Inset: Staggered magnetization m as a function of T at = 2. (b) 
Nonequilibrium phase diagram for a quench Ui ^ Uf with the fixed 
initial magnetization, m(0) = 0.4. For Uf > {<)Uf, the system fi- 
nally thermalizes to an AFM (PM) state. A nonthermal AFM order 
emerges in the colored region. The shading indicates the increasing 
life-time of the nonthermal AFM state as Ui is reduced. 



states in the single-site DMFT, the self-consistency condition 

'*'G<,(f,f') |2C 



reads K^(t,t') = t*^Go-(t,t') tMUM [A^(f,f'): hybridization 
function]. Since we are interested in the microscopic dynam- 
ics in a single magnetic domain, the system is assumed to take 
a spatially homogeneous configuration. 

In order to treat the long-time behavior of symmetry bro- 
ken states, we adopt the third-order weak-coupling expansion 
as an impurity solver, i.e., expand all the self-energy diagrams, 
including the Hartree term, by Weiss Green functions Goo-{t, t') 
(bare propagators) up to third order in U. Although the bare 
expansion is not a conserving approximation in the sense of 
Baym and KadanofF, it turns out to work remarkably well in 
the weak-coupling regime (U < 3). For instance, the total en- 
ergy is approximately conserved with negligibly small drifts. 
By comparison to quantum Monte Carlo (QMC) results lfl9ll . 
we confirmed that and m in equilibrium are correctly repro- 
duced [Fig.flla)], which is a considerable improvement from 
the Hartree approximation [Fig.lifa)] and the second-order it- 
erative perturbation theory II23I1 . 

Let us first look at results for quenches from f/, = 2 to var- 
ious Uf{< Ui). As shown in Fig.lJJa), m{t) quickly decreases 
after the quench due to the reduction of U, and starts to os- 
cillate coherently (amplitude mode) with a slow drift. As Uf 



decreases below ~ 1 .2, the oscillation disappears, and m ex- 
ponentially decays to zero. Assuming that the nonintegrable 
Hubbard model thermalizes, the long-time limit of the order 
parameter is determined by the thermal value mth at some ef- 
fective temperature T^ff. Since the total energy is conserved 
after the quench (f > tg) in the isolated system. Teg is given 
by the temperature of the equilibrium system with the same 
total energy. The final thermalized states are plotted as open 
dots in Fig.fTfa). Since we are considering rather slow changes 
(tg - 8) of U, the final states roughly keep track of the con- 
stant entropy curve [26J . 

The evaluated mth are indicated by aiTows in Fig.|2] and are 
plotted as a function of Uf in Fig. HI One notices that the cen- 
ter of the oscillation of m deviates more and more from mth as 
Uf is reduced. Surprisingly, at Uf = Uf - 1 .42, where mth 
vanishes cx — U]^f (Fig-IH with the mean-field exponent 
y6 = 5 (thermal phase transition), m still exhibits oscillations 
around a nonzero value for a long time. This suggests that 
the system is eff'ectively trapped in a nonequilibrium quasis- 
teady state, or close to a nonthermal fixed point, which al- 
lows for a long-Uved symmetry broken state with Teff above 
Tc- In the paramagnetic phase, the system shows prether- 
malization [9-11], i.e., momentum-integrated quantities such 
as the double occupancy thermaUze faster than momentum- 
dependent quantities (e. g., the momentum distribution). Here 
a new observation is that the order parameter m, even though it 
is momentum-integrated, also stays nonthermal, allowing the 
symmetry-broken state to survive for a long time. This can be 
attributed to the presence of 'classical fluctuations' [12] in the 
Hartree term, which is absent in the paramagnetic phase. 

To look at the qualitative change of the relaxation behav- 
ior around Uf ^ 1.2 more closely, we calculate the momen- 
tum distribution nt(r) = {cl.^^(t)cica-(t)) I I27I1 . In Fig. [3] one 
can clearly see the qualitative difference of «t between (a) 
Uf = 1.4 and (b) t/y = 1.2. In the former case, waves are 
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FIG. 3: (Color online). Time evolution of the momentum distribution 
rik for quenches (7, = 2 — >■ (a) Uf = 1.4 and (b) (// = 1.2. The red 
curves are thermal distributions achieved in the long-time limit. 



continuously generated at high energy, and cascade down to 
the lower energy region. They eventually reach the Fermi en- 
ergy ek - 0, and lead to an oscillation of the slope d^n at 
et = 12711 . In the latter case, the wave fronts never arrive at 
the Fermi energy but accumulate near = 0, which results in 
a steepening slope d^n. This evolution is opposite to a heat- 
ing effect, where an initially sharp momentum distribution is 
smeared out. Since the rtk in Fig. Ob) is very different from 
a thermal distribution (red curve) the fast relaxation of m for 
Uf < 1.2 [Fig.|2fa)] is due to dephasing, not thermalization. 

To characterize the nonthermal transition observed around 
{/ ~ 1 .2 quantitatively, we evaluate the relaxation time Tjeph 
for the dephasing of m(t) by fitting with e~'^'^'^^^. As shown 
in Fig. SJa), the dephasing critically slows down as Tdeph 
\Uf - t/;*|~' with t/"* = 1.23 (nonthermal transition point). 
At Uf = f/"*, m{t) shows a power-law decay of f^'''^ un- 
til thermalization starts to take place around t ~ 100. This 
evidences that one more quasi-critical point with an associ- 
ated diverging time scale exists away from the thermal crit- 
ical point {Uf - t/*). Moreover, a sharp kink is observed 
at Uf = t/"* in the plot of the inverse of the steepest slope 
(d,n)-^ = (max,{\dMt)\}r^ at ek ^ [Fig.gfa)]. Because 
a true discontinuity in the momentum distribution function, 
with (5(rn)^' - 0, would correspond to a power-law decay 
of the density correlations in space, one may thus note that 
at the nonthermal critical point the system evolves through 
an almost "critical state" before thermalization sets in. We 
also determined the frequency of the amplitude mode of 
m and the frequency (Ltg^,, of the oscillation of d^n at = 



for Uf > uf" by measuring the peak-to-dip distance of the 
oscillations. Note that near the critical point the period of the 
oscillation exceeds the life time (^ 100) of the trapped state, 
so that a meaningful measurement is not possible. However, 
the results in Fig.|4ja) indicate that w,„ and Luig^n extrapolate to 
zero as ~ \Uf — C/J* | . Based on this fact, we conclude that 
the amplitude mode is associated with the nonthermal fixed 
point, not with the thermal phase transition. This is not ex- 
pected in the Ginzburg-Landau picture, where the oscillation 
disappears when the curvature of the free energy potential at 
the origin changes sign at the thermal critical point. 

This quasi-critical point (or sharp crossover) becomes an 
exact critical point in the weak-correlation limit, where the 
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FIG. 4: (Color online). Various quantities used to characterize the 
qualitative change of the behavior around Uf = and Uf = (7"* 
(dashed lines) for quenches (/, Uf. Solid lines are guide for eye. 



dynamics is described by the Hartree approximation. As we 
show in Supplementary Material, the Hartree equation is, at 
half-filling, mathematically equivalent to the time-dependent 
BCS equation, which is known to be integrable with infinitely 
many conserved quantities [28, 29]. There is a strict transition 
for the motion of the order parameter from damped oscillation 
to overdamped decay that is both associated with a diverging 
dephasing time (overdamped decay) and a vanishing of the 
amplitude-mode frequency. What we found here is that the 
qualitative aspects of the transition are maintained even in the 
regime where the Hartree approximation breaks down {U > 
0.5) due to quantum corrections from higher-order diagrams. 
In fact, the Hartree equation gives quantitatively very different 
results in this regime [27]. 

As one increases f/,, the system spends less time near the 
nonthermal fixed point, and thermalization occurs earlier. For 
Ui - 2.5 [Fig.|2|b)], coherent amplitude oscillations are not 
visible anymore, and only a bump structure remains on a short 
time scale (f < 30) for Uf > 1.8. In this interaction regime 
the system does not show a clear signature of a transition, but 
a nonthermal crossover behavior is still seen in various quan- 
tities [Fig.Htb)] around Uf = f/"* ^ 1.85, which is estimated 
from the maximum of d^n. For Uf < t/"*, we find that the or- 
der parameter m shows a two-step relaxation [Fig.|5la)], i. e., 
the short-time and long-time dynamics have different expo- 
nential decay rates. The former is identified to be Tdeph, since 
it is smoothly connected to what we have defined as Tdeph in 
the previous f/, - 2 case. The latter is related to the thermal 
phase transition where mth disappears, hence denoted by Tth 
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FIG. 5: (Color online), (a) The log plot of m for quenches (/, = 
2.5 -s- Uf = 1.6, 1.7, 1.8, 1.9 from bottom to top. The black lines 
show fits of the two exponential relaxations, (b) The log-log plot of 
m for quenches Ui = 2.75 iff = 2.2, 2.25, . . . , 2.6 from bottom to 
top. The black line shows the slope of a power law decay cx 



density wave order if one translates the repulsive model to an 
attractive model [|32J. An open question of practical impor- 
tance is how to access this nonthermal fixed point. While we 
focused here on interation quenches, it will be interesting to 
consider also experimentally more realistic processes such as 
photo-doping induced by laser excitations. 

We thank H. Aoki, P. Barmettler, J. Berges, and T. Oka for 
fruitful discussions. The calculations were caiTied out on the 
Brutus cluster at ETH Zurich and on the UniFr cluster. We 
acknowledge support from the Swiss National Science Foun- 
dation (Grant PP0022-118866) and FP7/ERC starting grant 
No. 278023. 



SUPPLEMENTARY MATERIAL 



Hartree approximation 



POj. The obtained Tdeph and Tth are shown in Fig.|4tb). Inter- 
estingly, in most cases Tdeph is larger than Tth, that is, the slow 
dephasing of m is followed by faster thermalization. Further- 
more, thermalization is significantly delayed compared to Tth. 
At Uf = 1.9, for example, Tth = 15.3 while the delay time of 
thermalization is > 100. This allows the order parameter to 
survive longer than the thermalization time constant. 

Finally let us examine the relaxation around the thermal 
critical point. Thermalization critically slows down as one 
approaches the thermal critical point [Fig.HJb)] with 



TthOC |f//-f/*|-'> 



(1) 



which, unlike Tdeph, remains even when the interaction is in- 
creased. Since the critical behavior around the thermal tran- 
sition is universal, i.e., does not depend on details of the 
initial state or the the ramp protocol, it can be described 
by equilibrium properties. In fact, near a thermal (or quan- 
tum) critical point the relaxation time is known to behave as 
Tth ~ \Uf — fT3]. Here v is the critical exponent 

that characterizes the divergence of the correlation length, 
^ ^ \Uf — U^^\^^, and z is the dynamical critical expo- 
nent. Our result ([T]i is consistent with the mean-field expo- 
nents V = 5 and z - 2 for non-conserved order parameters 
lEEil]. Exactly at the thermal critical point (Uf = f/*), the 
order parameter thermalizes in a power law. In Fig.|5jb), we 
show the log-log plot of m around the thermal critical point 
(t/* - 2.40). The curve agrees very well with 

.-1/2 



m CX t 



(2) 



This can be interpreted as a sort of 'finite-size scaUng in time', 
i. e., an effect arising from the finite length of the real-time 
axis while the correlation time diverges. In fact, the relation 
(|2]l is consistent with m ^ t^fi/-^ with the mean-field exponent 
y6 = J, expected from the familiar finite-size scaling in space, 
m ^ L^l^l^ [L: linear system size) [31], with t U. 

We summarize our results in a nonequilibrium phase dia- 
gram in Fig. [lib)- Our results are applicable not only to an- 
tiferomagnetic order but also to superconductivity and charge 



In this section, we derive an equation of motion for the 
Hartree approximation, and show that it is mathematically 
eqtiivalent to the integrable time-dependent BCS equation 
ilH at half-filHng. 

Let us define a set of momentum distribution functions us- 
ing the nonequilibrium (lesser) Green function. 



nt{t)^{-i)Gt<{t,t) 



'ah</ 



ik-{R,-Rj) 



cl{t)cjAt)) 



(3) 
(4) 



with a,b - A,B sublattice indices and the number of sub- 
lattice sites. The nonequilibrium Green function satisfies the 
2x2 Dyson equation 



id, 

Sc 
6c 



^kcr ^kcr 

nBA nBB 



(5) 



and its conjugate equation 



r-AA 
r-BA 



Sc 
6c 



(6) 



where is the chemical potential, E° is the local self-energy 
on sublattice a, * denotes a convolution for time arguments, 
and 6c is the delta function defined on the Keldysh contour C. 
In the Hartree approximation, the self-energy is given by 



S^(f,f') = t/(f)4(0<5c(f,f'), 
^^Jt,t')^U(t)n^Jt)6c(t,t'), 



(7) 
(8) 



with n^(f) - {c]jt)cia-(t)) (i e a - A,B sublattice) the local 
density. In the presence of AFM order at half-filling, the local 
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densities are simply 



0.4 



n^it) = - + -o-m(t), 
n^t) = ^ - ^o-m(t). 



(9) 

(10) 



Using the Dyson equations Q and (|6]l with the Hartree ap- 
proximation, we obtain a closed set of equations of motion 
for the equal-time lesser Green functions, 

(id, + id,,)Gt^{t, f')|,'=, = [Gt^it, t) - Gil^t, t)] , (11) 

{id, + id,,)GZ.\t, t')\,,^, = -ek [Gff (f, t) - GZ^(t, t)] , 

(12) 

{id, + id,,)Gfi<{t, ?')!,'=, = [G^^(f, - G««<(f, 0] 

- U{t)m{t)crGf^<{t,t), (13) 

{id, + id,,)GZ<{t, f')|,'=, = -eu [G^"(?, - Gf^<(f, 0] 

+ U{f)m{f)&G^^l''{t,f). (14) 

Using Eq. (O to replace the equal-time Green functions by the 
corresponding momentum distribution functions, the equa- 
tions read 

d,[nt{t) + nZ{t)\^0, (15) 
d, [nZ{t) + ni^{t)] = -iU{t)m{t)cr [«^^(f) - ni^{t)] , (16) 
id, [«f^(0 - <(f)] = 2e, [nr^{t) - 4^{t)] 

+ U{t)m{t)cT[4i{t) + 4^{t)], (17) 
5, [n^(0 - <(0] = [nf^(0 - «^^(0] , (18) 

To make the expression transparent, we adopt a representation 
analogous to Anderson's pseudospin for the BCS theory Il33ll . 



(19) 



fl{t)^-Y,^i4a(t)~4a(t)i (20) 

/^(o = ^E^t«:^>)-<w]- (21) 

IT 

Then the Hartree equation can be simply written in the form 
of a 'Bloch equation', 

d,f,{t) = bk{t) X /,(f), (22) 

with a pseudospin 

fk-(fkJlfi) (23) 
and an effective magnetic field 

Mf) = (-2ei,0,f/(0m(f)). (24) 
The order parameter is self-consistently determined by 

m(f) = E/tXf). (25) 
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FIG. 6: Time evolution of m calculated with the Hartree approxima- 
tion for quenches (a) Uj = 2 — > Uf = 1.0, 1.1,..., 1.9 (from bottom 
to top) and (b) (/, = 2.5 Uf = 1.5, 1.6, . . . , 2.4 (from bottom to 
top). The color codes are the same as Fig. 2 in the main text. Note 
that the temperatures of the initial states are chosen to be different 
values from those of Fig. 2 in the main text to fix m(0) = 0.4. 



It turns out that the equation ( |22] ) is mathematically equivalent 
to the time-dependent BCS equation |28, 2^ if one appropri- 
ately translates the order parameter from the repulsive model 
to an attractive model |32]. The equation is known to be in- 
tegrable, and has infinitely many conserved quantities. For 



example, «j 
dospin' 



[Eq. (flST ll and the 'length of the pseu- 



lA I' ^(/i^)' + (/*)' + (/!) 



■z\2 



(26) 



are conserved for each k. 

Figure |6] shows the results obtained from the Hartree ap- 
proximation for the order parameter m. We choose the differ- 
ent temperatures of the initial states from those of Fig. 2 in 
the main text such that the initial value of the order parame- 
ter is the same as in the main text (m(0) = 0.4). For U, - 2 
[Fig.|6ja)], the overall tendency of the behavior is qualitatively 
similar to what we have observed with the nonequilibrium 
DMFT calculation in Fig. 2(a) of the main text, apart from 
the fact that the system never thermalizes within the Hartree 
calculation, such that the nonthermal transition becomes in- 
finitely sharp. After the quench, m rapidly decreases, and is 
suddenly trapped in a nonthermal value with coherent oscilla- 
tion of the amplitude. Between U j - 1.1 and 1.2, the center 
of the oscillations gradually approaches zero, and the behavior 
of m sharply changes from a damped oscillation to an expo- 
nential decay at some Uf - t/J*. However, quantitatively 
the evolution of m is very different. For example, the center 
of the oscillations, their frequency, the damping rate, and the 
transition point (t/"*) are all different from those of Fig. 2(a) 
in the main text. This is because the interaction strength that 
we consider here is already beyond the one {U < 0.5) for 
which the Hartree approximation works. It is thus surprising 
that the results including higher-order quantum corrections 
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FIG. 7: Time evolution of \ f^.\ at = calculated with nonequilib- 
rium DMFT for quenches (a) f/, = 2 -5> Uf = 1.0, 1.1, ... , 1.9 (from 
bottom to top) and (b) f/, = 2.5 Uf = 1.5, 1.6, . . . , 2.4 (from bot- 
tom to top). The color codes are the same as in Fig. 2 in the main 
text. The arrows indicate the corresponding thermal values reached 
in the long-time limit. 



shown in the main text nevertheless share the qualitative fea- 
tures of the Hartree approximation in this interaction regime. 
For Ui - 2.5 [Fig.|6jb)], the qualitative properties remain un- 
changed within the Hartree approximation, but sharpness of 
the transition is lost after quantum corrections are taken into 
account by nonequilibrium DMFT [Fig. 2(b) in the main text]. 

To see how the Hartree approximation fails for U > 0.5, 
we plot nonequilibrium DMFT results for |/j | (|26] i at et = 
(Fermi energy) in Fig.|7] These quantities would be conserved 
in the Hartree approximation. At e/c = 0, \f - \f^ \ since the 
off-diagonal Green functions (G^*, Gf^) are odd functions of 
ek- One can calculate /^|e=o from P"^ [Eq. ( |29] l below] us- 
ing the relation /~|,=o = 5 E<r(-'X^^ " K)- One can 



see in Fig.|7Ja) that is not conserved even for Uf - 1.0, 
but starts to decay immediately after the quench without any 
plateau. This suggests that the Hartree equation ( |22] | is not 
valid on any time scale, except in the very weakly correlated 
regime {U < 0.5). It was already clear from the equilibrium 
phase diagram [Fig. 1(a)], which shows the Hartree phase 
boundary as a dashed line, that there are large quantum cor- 
rections from higher order diagrams for U > 0.5. 



Momentum distribution function 

In this section, we show the derivation of the momentum 
distribution function 



rikit) 



-kiT 



{t)ckAt)), 



(27) 



and present the numerical results for the slope of the distribu- 
tion defik at e^^ = (Fermi energy), which sensitively measures 
whether and how thermalization takes place. 



By definition, one can obtain the momentum distribution 
(IZTT i from the nonequilibrium Green function (O, 



nka-{t) 



i^(-OGt^(M). 



(28) 



The Green function satisfies the 2x2 Dyson equation (|5]), 
which can be reduced to a set of 1 x 1 Dyson equations. To 
this end, we take the diagonal Green function at - 0, which 
we denote by 

n ^ GZU-o (29) 

(fl - A, B). It satisfies 

{i^, + ^i-K)*K-5c. (30) 

With f ^, the Green functions at arbitrary et are given by 



{i^, + ^i-K-4K)*GZ-5c, 



GAB _ ^ fAA J, n 



G 



ekO 



ka- 
kcr (T 



(31) 
(32) 
(33) 
(34) 



These are equivalent to solving a set of 1 x 1 Dyson equations. 

As is clear from Eqs. (I3ni-(l34li. the diagonal Green func- 
tions are even functions of et, while the off-diagonal Green 
functions are odd. Thus we have 



(35) 



To get the first derivative of the off-diagonal Green functions, 
we take a derivative with respect to eu and putting = in 
Eq. (|5]l to have 



{id,+^i-K)*d,Gt\, 
From this, we obtain 



■^ab I 
'ka-U 



=0 = 0. 



deGl%=0 



(36) 

(37) 
(38) 



As a result, the slope of the momentum distribution d^rtka- at 
ft = is calculated from a convolution, 

d,nUt)U=o ^-'-{P^^Pl + Pl* P^)"(f, f). (39) 

Numerical results for deti at = with the same param- 
eters as Fig. 2 in the main text are shown in Fig. [8] For 
Ui - 2.0 [Fig. |8ja)^ (b)], there are two clearly different be- 
haviors. When Uf > 1.2, deti coherently oscillates, which is 
because a wave mode created in the high energy region of the 
momentum distribution cascades down to the Fermi energy 
[Fig. 3(a) in the main text]. The slope d^n even changes its 
sign [Fig. Ob)] if the amplitude of the oscillations is strong 
enough. This never happens in thermal equilibrium. On 
longer time scales, the oscillation slowly damps, and the gra- 
dient d^n finally converges to the corresponding thermal value 
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FIG. 8; (Color online). Time evolution of the slope of the momentum 
distribution at = calculated with nonequilibrium DMFT 
for quenches (a) f/; = 2 f// = 1.0, 1.1, 1.2(< f/'*) (from top 
to bottom), (b) Ui = 2 ^ Uf = 1.3, . . . , 1.9(> f/^*) (from bottom 
to top), (c) Ui = 2.5 Uf = 1.5, 1.6, . . . , 1.8(< (/;*) (from top 
to bottom) and (d) Ui = 2.5 Uf = 1.9, . . . , 2.4(> Uf^) (from 
bottom to top). The color codes are the same as Fig. 2 in the main 
text. The arrows indicate the corresponding thermal values reached 
in the long-time limit. 



(thermalization). A sharp change occurs between Uf - 1.2 
and 1.3, where the amplitude of the oscillation in d^n is greatly 
enhanced, and even appears to diverge (x \Uf — with 
J/"* - 1 .23 [see the plot of the inverse of the steepest slope 
(5e«)~' = (maXf{|5en(f)|})^' in Fig. 4(a) of the main text]. 
As a result, a sharp jump of the momentum distribution starts 
to appear at the Fermi energy [Fig. 3(a) in the main text]. For 
Uf < 1.2, d^n overdamps without any oscillation. This qual- 
itative change of the behavior of dirn occurs at the same point 
(Uf - f/"*) as that of the order parameter m discussed in the 
main text. For f/, = 2.5 [Fig.[8fc), (d)], the enhancement of 
d^n is suppressed due to the limited life-time of the trapped 
state. However, one can still see a crossover of the relaxation 
behavior of from damped oscillation to overdamped decay 
around Uf ^ 1.85, the value of Uf for which the "steepest 
Fermi surface" is reached during the time evolution. Thus we 
use the steepest as a measure of the nonthermal transition 
point (f/"*) in the main text. 
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